% 定义Butcher表中的系数
A = [(88-7*sqrt(6))/360, (296-169*sqrt(6))/1800, (-2+3*sqrt(6))/225;
     (296+169*sqrt(6))/1800, (88+7*sqrt(6))/360, (-2-3*sqrt(6))/225;
     (16-sqrt(6))/36, (16+sqrt(6))/36, 1/9];
b = [(16-sqrt(6))/36; (16+sqrt(6))/36; 1/9];
c = [(4-sqrt(6))/10; (4+sqrt(6))/10; 1];

% 定义符号变量
syms z;

% 计算稳定性函数 R(z) = 1 + z * b' * inv(I - z*A) * ones(size(b))
I = eye(size(A));
ones_vec = ones(size(b));
R_z = 1 + z * b' * inv(I - z * A) * ones_vec;

% 化简表达式
R_z = simplify(R_z);
pretty(R_z)

% 计算极点（分母的根）
denominator = det(I - z*A);
poles = solve(denominator == 0, z);
disp('Poles of R(z):');
double(poles)

% 验证L-稳定性：计算 lim_{z->inf} R(z)
limit_R_inf = limit(R_z, z, inf);
disp(['Limit as z->inf: ', char(limit_R_inf)]);

% 绘制 |R(iy)| 在 y ∈ [-50, 50] 的图像验证A-稳定性
y = linspace(-50, 50, 1000);
R_iy = subs(R_z, z, 1i*y);
abs_R_iy = abs(double(R_iy));

figure;
plot(y, abs_R_iy);
xlabel('Imaginary axis y');
ylabel('|R(iy)|');
title('Stability function along imaginary axis');
grid on;
ylim([0 1.1]);
hold on;
plot(y, ones(size(y)), 'r--');  % 参考线 |R|=1